This invention is directed, in general, to simulating at least one characteristic of a physical system. In one aspect, the invention relates to a method for simulating a physical system such as a hydrocarbon-bearing reservoir to predict fluid properties and behavior in the reservoir.
Numerical simulation is widely used in industrial fields as a method of simulating a physical system by using a computer. In most cases, there is desire to model the transport processes occurring in the physical systems. What is being transported is typically mass, energy, momentum, or some combination thereof. By using numerical simulation, it is possible to reproduce and observe a physical phenomenon and to determine design parameters, without actual experiments using a model and apparatus. It can be expected therefore that design time and cost can be reduced considerably.
One type of simulation of great interest is a process of inferring the behavior of a real hydrocarbon-bearing reservoir from the performance of a model of that reservoir. The objective of reservoir simulation is to understand the complex chemical, physical, and fluid flow processes occurring in the reservoir sufficiently well to predict future behavior of the reservoir to maximize hydrocarbon recovery. Reservoir simulation often refers to the hydrodynamics of flow within a reservoir, but in a larger sense reservoir simulation can also refer to the total petroleum system which includes the reservoir, injection wells, production wells, surface flowlines, and surface processing facilities.
The principle of numerical simulation is to numerically solve equations describing a physical phenomenon by a computer. Such equations are generally ordinary differential equations and partial differential equations. As a means for numerically solving such equations, there are known the finite element method, the finite difference method, the finite volume method, and the like. In each of these methods, the physical system to be modeled is divided into smaller cells (a set of which is called a grid or mesh), and the state variables continuously changing in each cell are represented by sets of values for each cell. An original differential equation is replaced by a set of equations to express the fundamental principles of conservation of mass, energy, and/or momentum within each smaller unit or cells and of movement of mass, energy, and/or momentum between cells. These equations can number in the millions. Such replacement of continuously changing values by a finite number of values for each cell is called xe2x80x9cdiscretizationxe2x80x9d. In order to analyze a phenomenon changing in time, it is necessary to calculate physical quantities at discrete intervals of time called timesteps, irrespective of the continuously changing conditions as a function of time. Time-dependent modeling of the transport processes proceeds in a sequence of timesteps.
For most transport processes, the governing equations are nonlinear because the amount of mass, energy, or momentum in a cell and the movement of mass, energy, and momentum between cells typically have nonlinear relationships with the variables that define the physical state of the cells. In simulating a hydrocarbon reservoir, for example, the equations that model the reservoir are nonlinear partial differential equations that describe the unsteady-state flow of all fluids throughout the reservoir and relate the pressure and saturation changes of the fluids with time throughout the reservoir.
To simulate many physical systems, it is desirable to use implicit computations in which movement of a transported entity between cells depends on the solution at the end of a timestep. Implicit calculations require that the unknowns at the end of a timestep all be determined together. As a result, if the equations are nonlinear, the unknowns are typically computed using iteration. Iteration involves starting with some initial guess for the unknowns and applying some repetitive calculation to improve the guess until, after a sufficient number of iterations, the equations are satisfied to within some acceptable tolerance level. Since each iteration requires computing time, there is an incentive to use an iterative method that reduces the computing time as much as possible. Numerous iterative methods have been proposed for solving sets of nonlinear equations. One example is the well-known Newton-Raphson method.
The approximation used in a Newton-Raphson iteration results in a linear set of equations relating the unknowns at each cell to unknowns at its neighbors. These sets of equations are assembled into a global matrix equation that is then solved to obtain the next estimate of the solution. A similar matrix equation is obtained if the representation of the physical system is linear. In either case, the matrix equation is generally quite large and best solved iteratively. One iterative method for solving such matrix equations is a procedure called point Gauss-Seidel. In point Gauss-Seidel, a new solution estimate is calculated cell by cell. At each cell the new estimate is obtained by solving the mass, energy, and momentum balance equations for that cell, while holding unknown corresponding values at neighboring cells fixed at their latest estimates. In this procedure, a neighboring cell is one with which the current cell is in communication. A cell""s mass, energy, or momentum balance equations will contain terms multiplying the unknowns at its neighbors. When this calculation is repeated for all equations in the system, a new array of answers is created. This array is then checked to determine if the values satisfy the cell equations. To do this, it is convenient to define a residual (r) for each equation. If the new values satisfy the equations, then all residuals will be zero or very small. If not, the process is repeated with updated values of the unknown that are based on the previous iteration. This process is repeated until all the residuals are acceptably close to zero. This type of iterative method is called a xe2x80x9cpointxe2x80x9d iterative method because the method is performed a point, or cell, at a time.
Faster convergence can be obtained if point Gauss-Seidel is replaced by point-successive overrelaxation, or PSOR. In PSOR, the change in the estimated solution at each iteration is multiplied by an overrelaxation parameter, co, which must have a value between one and two.
Successful application of PSOR in simulation is generally limited to relatively simple models. Because PSOR methods are xe2x80x9cexplicitxe2x80x9d methods in which only one cell""s unknown values are calculated at a time, PSOR methods are prone to slow convergence. This shortcoming has led to efforts to include more implicitness in the solution methods. One method for doing this is called line-successive overrelaxation (LSOR). LSOR improves on PSOR by preserving implicitness in one direction. In it, the mass, energy, or momentum balance equations for a column or row of cells are solved simultaneously while the contributions of adjacent columns or rows are kept at their most recent estimates. Examples of LSOR applications can be found in (1) Mattax, C. C. and Dalton, R. L., Reservoir Simulation, Monograph Volume 13, Society of Petroleum Engineers (1990) and (2) Aziz, K. and Settari, A., Petroleum Reservoir Simulation, Applied Science Publishers Ltd, London, 1979.
The LSOR method used in the past has been applied primarily to models in which the cells are organized in a regular, structured grid having well defined rows or columns. Many models have been proposed with at least some of the cells arranged in a grid lacking this regular structure. It is believed that the practice of this invention represents the first application of LSOR principles to unstructured grids. Commercial use of unstructured grids has been slowed by the high cost of solution techniques for unstructured grids as compared to structured grids. A need exists for a simulation method that can be used to analyze representations of physical system using all types of cell configurations.
The method of this invention is used to simulate at least one characteristic of a physical system, regardless of whether the physical system has been discretized into cells occurring in structured or unstructured grids or a combination of both. The first step of the method is to discretized the physical system into a plurality of volumetric cells arranged adjacent to one another so as to have a boundary between each pair of neighboring cells. For each cell, linear equations are constructed that relate a cell""s state variables to the state variables of its neighboring cells. The next steps are to associate with each boundary a transportability value and then to rank the boundaries corresponding to a descending order of transportability values. The boundary rankings are then used to construct topologically one-dimensional strings of cells. A matrix equation is developed for each string by assembling the linear equations associated with the cells of each string. Improved estimates of the state variables of the cells are then obtained by solving the matrix equation for each string, one string at a time, until the matrix equations for all strings have been solved. This process is repeated iteratively until a convergence condition is satisfied. This solution produces state variables for all cells that simultaneously satisfy the linear equations for all cells. The state variables produced by the iteration can be used to simulate at least one characteristic of the physical system.
In a preferred embodiment, the construction of strings uses a rule that promotes the formation of strings having high transportability values at cell boundaries in the string cells.